
# empiricalMarginsNB(): This function is for calculating the change in expected counts and percentage changes in negative binomial regression models for various changes in a covariate of interest. It uses simulated coefficients from the estimated sampling distribution of the model input as "model_object" to calculate the 95% confidence intervals.

empiricalMarginsNB <- function(model_object, covariate, low, high, N){
  
  # Generate data frame allowing for fixed effects
  dat.mean <- dat.0 <- model.matrix(model_object)
  
  # Collect values
  x.min <- min(dat.0[,covariate])
  x.mean <- mean(dat.0[,covariate])
  x.median <- median(dat.0[,covariate])
  x.max <- max(dat.0[,covariate])
  x.sd <- sd(dat.0[,covariate])
  x.mean_sd_below <- max(c(x.mean - x.sd, 0))
  x.mean_sd_plus <- max(c(x.mean + x.sd, x.max))
  
  # Generate data frame matching values with names
  values <- c(x.min, x.mean, x.median, x.max, x.mean_sd_below, x.mean_sd_plus)
  names <- c("min","mean","median","max","mean_sd_below","mean_sd_plus")
  
  # Put in custom values of covariate
  dat.0[,covariate] <- v1 <- values[names == low]
  dat.mean[,covariate] <- v2 <- values[names == high]
  
  # Simulate beta's from estimated sampling distribution of model, transpose
  beta_dist <- as.matrix(t(mvrnorm(N,
                                   mu = coef(model_object),
                                   Sigma = vcov(model_object))))
  
  # Calculate predicted values and difference
  diffs <- c()
  pcts <- c()
  for(i in 1:ncol(beta_dist)){
    p.count.1 <- exp(dat.0%*%beta_dist[,i])
    p.count.2 <- exp(dat.mean%*%beta_dist[,i])
    diffs[i] <- mean(p.count.2 - p.count.1, na.rm = T)
    pcts[i]<- mean(((p.count.2-p.count.1)/p.count.1)*100, na.rm = T)
  }
  
  # Return mean marginal effect in sample with 95% confidence interval
  return(rbind(c(effect = mean(diffs),
                 quantile(diffs, c(0.025, 0.975))),
               c(pct_change = mean(pcts),
                 quantile(pcts, c(0.025, 0.975)))))
}






# twostageboot(): This function estimates coefficients and bootstrapped confidence intervals of endogenous regressors in two-stage models where the first stage coefficients are estimated via OLS. The function requires the two models to be estimated separately prior to the use of the twostageboot() function. 

twostageboot <- function(B, model_1, model_2, regressor, base_data, log_t = F){
  require(foreach)
  require(dplyr)
  
  ## Loop to extract second-stage coefficients
  coefs <- foreach(i = 1:B, .combine = "c") %do% {
    b_data <- sample_n(base_data, 
                       size = nrow(base_data), 
                       replace = TRUE)
    
    ## Stage 1
    stage_1 <- update(model_1, 
                      data = b_data)
    b_data[,regressor] <- predict(stage_1)
    if(log_t == T) {b_data[,regressor] <- (1/(1+exp(-predict(stage_1))))*100}

    
    ## Stage 2
    stage_2 <- update(model_2, data = b_data)
    coef(stage_2)[regressor]
  }
  
  se <- sd(coefs)
  return(c(beta = mean(coefs), 
           quantile(coefs, c(0.05, 0.95))))
}





